On the role of metastable intermediate states in the homogeneous nucleation of solids 

from solution^! 
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The role of metastable liquid phases in vapor-crystal nucleation is studied using Density Functional 
Theory (DFT). The bulk free energy functional is modeled using a hybrid model consisting of a sum of 
a hard-sphere functional (Fundamental Measure Theory) and a liquid-like contribution to account for 
the attractive part of the interactions. The model gives a semi-quantitatively accurate description of 
both the vapor-liquid-solid phase diagram for both simple fluids (Lennard- Jones interactions) and of 
the low-density/high-density/crystal phase diagram for model globular proteins (ten Wolde-Frenkel 
interaction) . The density profile is characterized by two local order parameters, the average density 
and the crystallinity. The bulk free energy model is supplemented by squared-derivative terms in 
these order parameters to account for inhomogeneities thus producing a model similar in spirit to 
phase-field theory. It is shown that for both interaction models, the vapor-crystal part of the phase- 
diagram can be separated into regions for which metastable liquid phases are more or less stable 
than the vapor, but always less stable than the solid. The former case allows for the possibility of 
double nucleation whereby liquid droplets nucleate from the vapor followed by a separate nucleation 
of the solid phase within the liquid droplets. Whether or not this actually occurs depends on the 
relative free energy barriers for vapor-solid and vapor-liquid nucleation and it is shown that for 
simple fluids, double nucleation is indeed favored at sufficiently large supersaturation. Finally, by 
studying the minimum free energy path from the vapor to the solid, the separate possibility of 
transient nucleation, where the vapor-solid transition involves a single nucleation event but the sub- 
critical clusters tend to be liquid-like while small, is shown to also be possible when the metastable 
liquid exists, but the supersaturation is too small for double nucleation. 
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I. INTRODUCTION 

The goal of this work is to describe, at a semi-microscopic level, the process of homogeneous nucleation of solids 
from solution. This is studied in the approximation that the effects of the solvent can be accounted for by an effective 
interaction between the solvate molecules so that at low concentrations, the molecules in solution are treated as a 
single-component gas with the crystalline phase being the corresponding solid phase in the effective single component 
system. In this picture, there is also the possibility of a dense liquid-like phase which simply corresponds to the 
liquid phase of the effective single-component system. Thus, the following work is intended to be applicable both to 
nucleation in single-component systems, such as simple fluids, and to nucleation of a solid phase from macromolecules 
in solution. For the sake of consistency, most of the discussion will be phrased in the language of nulceation of solids in 
a single component system. So, in the following, all statements concerning the "vapor phase" should be understood to 
be equally applicable to the "low-concentration" or "low density" phase of molecules, such as proteins, in solution and 
statements concering the "liquid phase" are equally applicable to the "high-concentration" or "high-density" phase 
of a solution. 

In general, homogeneous nucleation involves at least two phases: the stable phase, S, and the metastable phase, M. 
Initially, the system consists of M in its homogeneous, bulk state. Fluctuations give rise to small volumes of S which 
will here always be taken to be spherical clusters. By definition, the free energy of bulk S is lower than that of bulk 
M but finite clusters also involve an interface which has higher free energy, on average, than either phase. Hence, 
sufficiently small clusters are typically unstable, having higher free energy even than the equivalent amount of M. 
Furthermore, for small clusters, there is no reason to believe that the material inside the cluster is in the bulk S state 
since the interface can have a volume comparable to that of the bulk region. Thus, small clusters will be unstable 
with respect to M: indeed, will be unstable with respect to smaller clusters thus leading to cluster dissolution. For 
sufficiently large clusters, the interior will consist mostly of material near the bulk S state, and the volume of the 
interface will be small compared to the volume of the bulk so that the cluster is stable with respect to M but unstable 
with respect to cluster growth thus driving the transition of the whole system from M to S (see, e.g. 1|). 
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itn respect to cluster growth thus driving the transition ot the whole system from M to b (see, e.g. 
In recent years, it has become apparent, first for proteins0-@ but then even for simple fluids]!, 0, 
geneous nucleation of solids from vapor can be more complicated than this simple picture. This is because in many 
cases, the conditions for solid nucleation also allow for a second metastable phase, namely that of the bulk liquid. In 
cases where the bulk liquid is more stable than the bulk vapor (i.e. has lower free energy) but less stable than the 
solid, it can be that, instead of directly nucleating a crystalline cluster from the vapor, it is energetically favorable 
to first nucleate a liquid cluster which then grows and, subsequently, to complete the transition to the crystalline 
state via a second homogeneous (liquid-crystal) nucleation event within the liquid cluster (thus creating a growing 
solid cluster within the growing liquid cluster). This process will be referred to as "double nucleation". Whether 
or not this actually occurs will depend on the free energy barrier encountered in the vapor-crystal and vapor-liquid 
transition. Alternatively, even if direct vapor-crystal nucleation is favored, it is still possible that the growing cluster 
nevertheless passes through a liquid-like stage which, however, is always sub-critical, as a transient on the way to 
creating a solid-like critical cluster. In this case, small clusters would be liquid-like in structure, becoming more 
crystalline as they grow larger. This scenario, involving only a single nucleation event, will be referred to as "tran- 
sient two-step nucleation" . The simulations of ten Wolde and FrenkelQ are clear illustrations of transient two-step 
nucleation: they exhibit a nucleation pathway that involves liquid-like small clusters followed by solidification as the 
cluster approaches the critical size. Vekilov discusses both scenarios and suggests that even when double nucleation 
is possible and energetically favored, it may be suppressed by kinetics^]. Kaschiev and Vekilov have analyzed the 
effect of double nucleation on observed nucleation rates 0. van Meel et al report simulation results showing double 
nucleation for a Lennard- Jones fluid Q. 
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The key to understanding these processes is the construction of models for the free energy of inhomogeneous - 
multiphase- systems. Indeed, Classical nucleation theory (CNT), from which we take our lead, is fundamentally 
based on a description of the free energy as a function of the size of the cluster. In CNT (which should be thought of 
in terms of the simpler vapor-liquid transition) the width of the cluster interface is taken to be zero and the interior 
of the cluster is assumed to be in the bulk state, S. The only variable is then the size of the cluster and it is assumed 
that homogeneous nucleation consists of the growth of a cluster from size zero until it is arbitrarily largefTol fTl| . This, 
in other words, is understood to be the "nucleation pathway". This idea can be developed more rigorously to include 
finite width of the interface and non-bulk properties of the interior in which case the nucleation pathway involves a 
description of the variation of all these properties - size, interfacial width and interior density - simultaneously fl2| . 

The process of solid nucleation is more complex, as the phenomenology above would imply. First, the vapor-solid 
transition involves at least two order parameters: the density and the crvstallinitv [H. IH [l3l Il4| . Since the typical solid 
density is close to that of the metastable liquid (if it exists), the double nucleation scenario involves a separation of 
these order parameters: the vapor-liquid transition involves two different densities but occurs at zero crystallinity while 
the liquid-solid transition occurs at (nearly constant) density and involves a change from zero to finite crystallinity. 
In the transient two-step scenario, a small (unstable) cluster begins to form at zero crystallinity but at some point, 
as the cluster becomes larger, the crystallinity increases so that the critical cluster has both finite density difference 
compared to the vapor and finite crystallinity. Both of these can be viewed as "two-step" scenarios compared to the 
possibility that, from the beginning, the crystallinity and density change at the same time. Clearly, description of 
these processes must involve a model for the free energy of inhomogeneous systems which is sufficiently detailed so as 
to describe liquid and solid clusters embedded in a vapor bulk. In order to capture the transient scenario, this must 
be supplemented by some means of determining the nucleation pathway. 

Here, the basic tool for calculating the free energies will be classical, finite-temperature Density Functional Theory 
(DFT). If DFT were sufficiently well-developed and technically simple, nothing more would need to be said about this 
part of the study: unfortunately, neither of these conditions holds for the solid phase. (In contrast, direct calculations 
for liquid- vapor systems are possible [ill. Il5j . DFT is only sufficiently well understood so as to give an a priori 
description of bulk fluid, bulk solid and inhomogeneous fluid-solid systems for hard-sphere interactions. For any more 
realistic potential - particularly those with attractive interactions - some sort of modelling is necessary. Typically, this 
will involve the introduction of an effective hard-sphere diameter and the representation of the free energy as the sum 
of the hard-sphere functional and a second term accounting for the attractive part of the interactions. Futhcrrmorc, 
even for the hard-sphere functional, the calculation of the free energy for an inhomogeneous system (i.e. of a solid 
cluster embedded in a fluid background) is very complicated and computationally expensive [la]. Fortunately, an 
easier alternative is available. It has been shown that the exact free energy for a solid can be expanded in terms of 
gradients of the order parameters [l7l - [l9j thus providing a connection between DFT and the older gradient theories of 
Cahn et al(2^,[2l|, as well as the phase-field theories commonly used to study solid-solid transitions [22j|. This is the 
justification for using a gradient model in the present work. A significant advantage is that the free energy functional 
is then only needed for the case of homogeneous (i.e. bulk) systems, thus placing less stress on the accuracry of the 
DFT model. On the other hand, expressions for the coefficients of the gradient terms must be found. In principle, the 
exact results express these in terms of derivatives of the bulk free energy but in practice, they are hard to calculate 
except for the case of a fluid. In this work, a semi-empirical procedure will be used to fix these terms. The main 
difference between the present work and that of Gunton et alQ, which is similar in spirit, is that the latter made use 
of toy free energy models whereas here realistic models are used which give semi-quantitatively accurate bulk phase 
diagrams as well as of the liquid- vapor phase transition[l5j. 

Given these approximations, the free energy for any configuration of order parameters can be calculated. The 
practical exploration of the models will make use of energy-surface techniques commonly applied to the study of 
chemical reaction pathways and structural transitions (23|. A recently developed method - involving the approximation 
of the spatially-varying order parameters as piecewise-continuous functions, will be used to determine the critical 
clusters - i.e. saddle points of the free energy - for homogeneous nucleation. This will already give enough information 
to identify if and when double nucleation is possible. The nucleation pathway will be identified, as is commonly done 
for chemical reaction pathways, using steepest descent paths. These identify the most likely path for the transition 
given the free energy surface and are a natural generalization of the CNT pathway. They do not take into account 
kinetic factors, such as rates of mass transport, that could play a significant role particularly for small molecules. 

Section II will describe the technical details of the free energy models used. The bulk thermodynamics is used in 
Section III to limit the regions of the phase diagram in which double nucleation is possible. A simple model for double 
nucleation is also used to illustrate the role of bulk free energy differences and of surface tension. The fourth Section 
describes the application of the model to planar interfaces and illustrates the role of wetting. Detailed calculations 
of the energy barriers for direct nucleation of the crystal and of double nucleation are presented in Section V. The 
possibility of transient double nucleation is also described. The paper ends with a summary and with a discussion of 
future directions. 
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II. THE FREE ENERGY MODEL 
A. Density 

The central quantity denning the system state in DFT is the local density, p(r). It is possible to formulate the 
theory with no a priori restrictions placed on the density, but this is computationally expensive. A commonly used 
alternative is to represent the density in terms of a set of basis functions. For bulk crystalline systems, this usually 
means a sum of Gaussians localized at the lattice sites, 

p(r;x, piatt, a) = (—^ exp (-a (r - R 4 (pi a tt)) 2 ^j (1) 

i 

where Ri is a lattice vector, a controls the width, x is the occupancy and pi a tt is the lattice density. The average 
density is 

P=yJ P (r;S, Piatt, a) = x pi att (2) 
The density can also be written in Fourier representation as 

p (r; x, piatt, a) = xp ia tt ^ exp (ir ■ K, (pi att )) exp (-K 4 2 (piatt) /4a) (3) 

i 

where K, is a reciprocal lattice vector. Notice that in this representation, it is clear that the limit a — > gives a 
uniform density p(r) = ~p which describes a fluid. It is therefore natural to characterize the crystallinity by the size of 
the amplitude of a typical nonzero wavevector term such as 

X = exp(-K? (pi a tt)l±a) 

so that the density can also be written as 

p (r; x, p latu a) = xpi at t J2 ex P (* r ' K * (Piatt)) x (^(«»«)/^(«»t t )) 2 (4) 

i=0 

thus showing that the density is parameterized entirely by x, piatt and \- I n the following, we neglect the variation 
of the lattice density and generalize to inhomogeneous systems by allowing the occupancy and the crystallinity to 
depend on position. It is also convenient then to replace the occupancy by the average density so that we have 

p(v) =p(r)^exp(zr-K l (p/ att ))(x(r)) (Kl(pi " tt)/Kl(p, "" ))2 (5) 

i=0 

The order parameters are then the local average density, p (r) , and the crystallinity, \ ( r ) • It will sometimes be more 
convenient to replace the latter by the amplitude of the smallest non-zero wavevector, pi(r) = p(r) \ (r). 

B. Gradient expansion 

In order to determine the density, a model for the (grand canonical) free energy functional, il[p] is necessary. Good 
models exist for liquids and can be used to study, e.g., the liquid- vapor transition in great detail[l5j. However, 
the theory is less well developed for the solid phase and in any case calculations for inhomogeneous solids are very 
expensive. The present work therefore makes use of a gradient expansion of the free energy which focusses attention 
on the order parameters and only requires information about homogeneous solids [TlHll-" 

The grand-canonical free 

energy is written as 

fl[p] = F[p] -pj p(v)dv (6) 

where p is the chemical potential. In general, if the density can be written in terms of n order parameters, T = 
{T 1: ...,r„} , as 



p(r) = /(r;r(r)) 



(7) 
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so that the density of the uniform, bulk system is 

pr(r)=/(r;r) (8) 
and if T (r) in some sense "slowly varying", then the squared-gradient approximation (SGA) for the functional F[p] is 

FH ./{/(r W > + I™>^^} ( 9 , 

where the summation convention is used. The first term of the free energy involves the bulk free energy density 
defined as 

/(D = lF[ Pr ] (10) 
For the order parameters used here, the free energy is explicitly 



F[p] 



f (p (r) , X (r)) + \KZ (p (r) , X (r)) Mil Mil 1 
+K% (p (r) )X (r))|?¥ + |^ 6 X (p (r),x(r))^^ ' 



C. Bulk Free Energy 

The bulk free energy model used here is based on the idea of separating the free energy into a hard-sphere contri- 
bution, for which the DFT is well developed, and a second contribution that accounts for the long-ranged attractive 
interactions. A particularly simple model is based on the observation that the local structure of an FCC solid and 
a simple fluid are quite similar so that, as a simplest approximation, one can imagine that the correction to the 
hard-sphere model is independent of the local structurejj, [24| giving 

±F [ Pr ] = 1f hs [ Pr d (r) 3 ] + f tail (p r> d IT)) (12) 

where d IT) = d (p) is the effective hard-sphere diameter and 

ft a a (Pr,d (T)) = ±F [p] - ±F HS [pd 3 ] (13) 

is the contribution of the attractive part of the interaction to the liquid free energy. The effective hard-sphere diameter 
is determine using the WCA expression as modified by Ree et alj25l. [26j. 

D. Bulk phase diagram 

In DFT, one is always working in the grand-canonical ensemble so the external parameters are the temperature, 
the chemical potential, p and the applied external field, 0(r). The latter includes any confining walls: if the walls 
are hard, then the volume is a fixed parameter (as will always be assumed here). The appropriate free energy is the 
grand potential, 



tt\p] = F\p]-pL J p(r)dr+ J 0(r)p(r) dr 



and it should be noted that all information about the state is encoded in the local density function, p (r) . The 
equilibrium states (i.e. density distributions) are determined by minimization, 

0=™ (») 
dp (r) 

which is to say 



(15) 
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For parameterized profiles, the requirement that the free energy be a minimum gives 



and if the parameters are constants, then 



n Sn ^ Ha 



o = ^ (i 7 ) 



For a given value of chemical potential, there may be multiple solutions for the density in which case the equilibrium 
state is the one corresponding to the absolute minimum of the grand potential. Two phase coexistence therefore 
requires that 



SF[p] 



H-(j)(r)- 



6 P (r) 

n [pi] = n [p 2 ] 

In particular, using the chain rule for functional differentiation, 

dF [p] f 6F [p] dp (r) 



(18) 

pi 



dp J Sp (r) dp 
6F[p] 



dr (19) 



Sp (r) 



dr 



gives the usual thermodynamic relation 



dp ^ V 

Thus, for uniform densities, pi (r) = 'p i and e.g. F [p^ — > F (Pi), the conditions for coexistence are 



d v F[p] p-1 / 0(r) rfr . 



3- Px »-vl^ d *=^wr (21) 

n(pi) = n(pa) 

which is to say equality of chemical potentials and of pressures since in the bulk, Q = —PV. 

For the crystalline system, the parameterization used here involves not just the average density but also the 
crystallinity and the lattice parameter. The conditions for an extremum of the free energy are then 



dyF (p,x, Piatt) 

dp 


1 d 

-'^v&p 1 


dy F i~P,X,Piatt) 


i d r 


dx 


~vdx~J 


9yF(p,x, Piatt) 


1 d 


dpiatt 


V dpiatt 



J p(r;p, X , Piatt) <f>(r)dr (22) 



p{r;p,X, Piatt) 0(r) dr 

In principle, the chemical potential (an external held) arc specified and these equations solved for p,x and piatt- in 
practice, it is simpler to choose a value of the lattice density and to use these conditions to determine p, x and p. 
Further technical details are discussed in Appendix [X] 

III. THERMODYNAMICS OF TWO STEP NUCLEATION 
A. Independent variables and Ensembles 



The main calculational tool used here, DFT, is formulated in the grand-canonical ensemble in which the independent 
variables are the temperature, chemical potential and volume (or, more generally, applied held) and the free energy 
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of interest is the grand potential. Experiments are typically performed at constant pressure, temperature and volume 
for which the Gibbs free energy is relevant. Fortunately, in discussing nucleation, one is always interested in free 
energy differences and these are independent of the ensemble]!?]]. In the following, since it is based on DFT, the 
grand-canonical ensemble is always used, however in many cases the difference between ensembles can be suppressed 
by focussing on physical quantities. Thus, rather than specify the chemical potential (for the grand ensemble) or the 
pressure (for the PVT ensemble) one can specify the density of the final phase which is a meaningful variable in both 
formulations and in either case is uniquely related to the state variable. 



B. Interaction potentials and the fluid phase 



In this study, the nucleation properties of two different systems will be considered: simple fluids as described by 
the Lennard- Jones potential, 

•~M -*((?)"- (7)0 < 23 » 

and globular proteins as described by the ten Wolde-Frenkel model interaction, 

oo, r > a 

vtw F (r) = { 4e ( ( i \ 6 „ ( i ^ ... , (24) 



which will be studied here for the value a = 50 as is appropriate to model the phase behavior of globular proteins. 

Both of these potentials are long-ranged in the sense of decaying as power-laws. In simulation, infinite-ranged 
potentials arc difficult to deal with, so any long-ranged potential v(r) is typically cutoff at some distance, r c . For 
Monte-Carlo, a simple shift of the potential to make the energy continuous at the cutoff is typically used so that the 
so-called truncated and shifted potential is 

, . / v(r) - v (r c ) , r < r c . > 

«mc(r) = | 0r>rc (25) 

In molecular dynamics simulations, the force is usually required to be continuous so that the force-shifted potential 
is commonly used, 

,, m _ / v ( r ) ~ v ( r c) ~ v' (r c ) (r-r c ), r < r c , , 

Vmd v>-\ 0, r > r c [M> 

where v'(r) = dv(r)/dr. Other forms are also important, particularly the Broughton-Gilmer modification of the LJ 
potential: 

VLj(r)+d, r<2.3a 

vbg (r) = { C 2 (f ) 12 + C 3 (f ) 6 + C 4 (f f + C B , 2.3a < r < 2.5a (27) 

0, r > 2.5a 

where d = 0.0I6I32e, C 2 = 313. 66e, C 3 = -68.069e, C 4 = 0.083312e and C 5 = 0.74689e[28j. 

All of these details significantly affect the liquid-vapor phase diagram. The DFT model described above requires 
as input both the interaction potential and the equation of state for the fluid phase. For the LJ potential with no 
cutoff, essentially exact empirical equations of state are available for temperatures above the triple point [29l - [3l| . For 
finite cutoffs, these must be modified with inexact mean-field corrections. Although useful for studying liquid-vapor 
coexistence [ill EH; they are of limited utility for vapor-solid nucleation since the interesting affects occur below the 
triple point. Less accurate, but more broadly applicable, are mean-field equations of state based on thermodynamic 
perturbation theory such as the Weeks-Chandler- Andersen (WCA) theoryJH, [2(| l32l - l3"H which will be used here. 



C. Solid phase diagrams and intermediate phases 



Figure [T] shows the calculated vapor-liquid-FCC solid phase diagram for the infinite-ranged LJ potential and for 
the truncated and shifted tWF potential compared to simulation in terms of the reduced temperature T* = ksT/e. 
The LJ phase diagram possesses both a liquid-vapor critical point and a triple point whereas for the protein model, 
the liquid-vapor transition is metastable. Since the fluid equation of state is being calculated from thermodynamic 
perturbation theory, which is a mean field theory, the critical point is, as expected, poorly described for both potentials. 
Apart from this expected inaccuracy, the model is in semi-quantitative agreement with simulation for all three phases. 
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1. Nucleation scenarios for globular proteins 



In order to clarify the thermodynamics of mctastable states, we consider in more detail the phase diagram derived 
from the tWF potential. Figure [2] shows a line crossing the phase diagram at constant temperature. The points A 
and B identify coexisting vapor and solid phases which, by definition, have the same free energy and the same control 
parameter (chemical potential or pressure). Starting at the vapor point A and moving along the isotherm in the 
direction of increasing density, i.e. to the right, corresponds to increasing the chemical potential and to each value 
of the density there will be a point on the isotherm to the right of the solid-branch of the binodal, i.e. to the right 
of point B, having the same chemical potential; however, under these conditions the solid phase will have lower free 
energy than the vapor phase so that the vapor points to the right of A are metastable with respect to the solid phase. 
Eventually, as one moves along the isotherm, the point in the vapor region reaches the vapor-liquid coexistence curve 




25 1 1 1 1 1 

0.5 1 

FIG. 2. The phase diagram for the tWF potential as in Fig. [T]but without the simulation data and showing the spinodal for the 
vapor-liquid transition (broken line). The figure also shows, in red, the boundaries for the existence of the intermediate liquid 
phase: there is no liquid phase for chemical potentials corresponding to vapor densities to the left of the red line. A similar 
line for the solid phase is nearly indistinguishable from the solid binodal. Similar boundaries also exist on the hennard- Jones 
phase diagram but are so close to the binodals as to be almost indistinguishable. The horizontal line, an isotherm, picks out 
coexisting states with the coexisting vapor and solid states being labeled A and B, respectively. 
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so that there is a coexisting liquid phase. However, by definition, this particular liquid state has the same free energy 
as the vapor, one knows that it has higher free energy than the solid and so can only be metastable. The remains true 
as one moves along the isotherm to densities to the right of the vapor branch of the liquid- vapor binodal. Eventually, 
one reaches the vapor-liquid spinodal and above this density, the vapor phase does not exist. The set of liquid points 
with the same chemical potential as the vapor points on the spinodal will therefore divide region to the right of the 
liquid binodal into a lower-density region, for which a vapor with the same chemical potential can always be found, 
and a higher-density region with no corresponding vapor. 

What is the role of the metastable liquid phase in vapor-solid nuclcation? We can only address this question here 
with respect to double nuclcation. (Transient metastable states are a property of the nuclcation pathway and its 
presence or absence cannot be answered solely from a consideration of the bulk phases.) First, consider again the 
vapor-liquid coexistence curve. By a reversal of the reasoning above, vapor states to the left of the vapor branch 
correspond to liquid points to the left of the liquid-branch of the coexistence curve and are more stable than the 
corresponding liquid state. Moving to the left away from the vapor branch, the corresponding liquid point also moves 
left until it reaches the liquid-vapor spinodal: beyond this point, there is no corresponding liquid. This therefore 
defines a line in the phase diagram with the property that vapor states to the left have no corresponding liquid states, 
shown as a red line in Fig. [5J In particular, points on the vapor branch of the vapor-solid coexistence curve have 
no corresponding liquid. This line therefore divides the metastable vapor region into two parts: a low density part 
in which there is no corresponding liquid state (i.e. no liquid with the same chemical potential) that can play a role 
in nucleation and a higher density region where the corresponding liquid exists as a metastable state. Hence, for 
systems prepared with the vapor density between the vapor-branch of the solid-vapor coexistence curve and the new 
demarcation line, double nucleation is not possible as there is no metastable liquid. The region within which vapor 
and liquid states with the same chemical potential can be found is therefore an envelope around the binodal and will 
be referred to as the /x-nodal curve. 

Starting at the vapor branch of the vapor-solid coexistence curves, and moving to the right we therefore find that: 

1. From coexistence to the /x-nodal line, there is no liquid phase and so, double nucleation is not possible. In this 
region, the solid free energy is less than that of the vapor. 

2. From the /x-nodal line, to the vapor branch of the vapor-liquid coexistence curves, there is a liquid state but 
it has higher free energy than the vapor (which has higher free energy than the solid). Nucleation via a liquid 
cluster is possible, with the solid cluster nucleating within the liquid but the liquid cluster would always be 
metastable. This is a form of transient two-step nucleation. 

3. From the vapor branch of the liquid-vapor coexistence curve to the vapor-liquid spinodal, the liquid has lower 
free energy than the vapor but higher than the solid. True double nucleation is possible depending on the 
barriers for the formation of a critical liquid cluster within the vapor, a critical solid cluster in the liquid as 
compared to the barrier for directly forming a critical solid cluster in the vapor. Even if the latter is favored, a 
transient scenario is possible. 

All of this is shown in a more direct way for the tWF potential in Fig. [3] where supersaturation, S = P v /P C oex 
with P v the vapor pressure and P CO ex the vapor-solid coexistence pressure, is used as the independent variable. By 
definition, for S < 1, the vapor is the stable phase and for S > 1 the solid is the stable phase. The latter region 
is divided into three sections: that for which there is no corresponding liquid state, that in which there is a liquid 
but it has higher free energy than the vapor and that in which there is a liquid state with lower free energy than the 
vapor but higher free energy than the solid. Only in the latter case is liquid nucleation possible so it is only in this 
region of the phase diagram that double nucleation can occur. This therefore represents a set of necessary conditions 
for double nucleation. Nucleation pathways that do not involve double nucleation but which pass through liquid-like 
state, i.e. transient two-step nucleation, can in principle occur for any value of S > 1. The primary goal of detailed 
DFT calculations is to determine the necessary conditions for double nucleation and to assess for what conditions, if 
any, transient two-step nucleation occurs. 

2. Nucleation scenarios for simple fluids 

In some ways, the phase diagram for simple liquids is more complex because the various coexistence curves cross. 
As shown in FigJ5J the calculations indicate that there is liquid-vapor coexistence below the triple point. However, 
just as in the case of the globular proteins, the vapor branch of the vapor-liquid coexistence curve is to the right 
of the vapor branch of the vapor-solid coexistence curves thus indicating that the coexisting liquid and vapor have 
higher free energy than does the solid phase (at the value of chemical potential or pressure corresponding to liquid- 
vapor coexistence). Hence, below the triple point, the liquid phase is again metastable just as in the case of the 
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FIG. 4. Phase diagram for the LJ potential with (vapor-solid) supersaturation as the independent variable. 



globular proteins. Similarly, below the triple point the liquid branch of the liquid-solid coexistence curves lies to 
the left of the liquid branch of the liquid-vapor coexistence curve: i.e. it is in the metastable (or even unstable) 
region of the liquid-vapor phase diagram. So, these liquid states have higher free energy than the corresponding 
vapor phase and the coexisting liquid and solid phases are metastable with respect to the vapor. Taking all of this 
together, only the vapor-solid transition is usually drawn below the triple point, but there is a metastable liquid and 
a liquid-vapor coexistence curve in this region. It is this metastable liquid phase which gives rise to the possibility of 
two-step nucleation even for simple liquids. However, the fact that the vapor branches of the vapor-solid coexistence 
curves and the vapor-liquid coexistence curves are very close together makes it more difficult to display the various 
metastability boundaries. Figure [4] shows the phase information in terms of the supersaturation and the similarities 
and differences to the model protein behavior are evident. The main difference is that the intermediate liquid state 
exists for all values of the supersaturation leaving only the two regions distinguishing liquids which are less or more 
stable than the vapor phase. 
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D. A simple picture of double nucleation 

In this subsection, the goal is to anticipate the following, more technical, developments to give an idea of how 
double-nucleation can be described by something that looks like an extension of CNT. We begin with the free energy 
functional defined above specialized to spherical symmetry, 

OH-/ {« (?M.*MhHW + *~^^ + ^W} (28) 

where u = / — pp is the grand potential per unit volume and the coefficients of the gradient terms are taken to be 
constants. It is assumed that the temperature and chemical potential are such that there are three bulk states: a 
solid with order parameters (p s ,Xs)j a liquid with order parameters (p;,0) and a vapor with parameters (p, u ,0).We 
now introduce a very simple model for the order parameters for a (spherically symmetric) vapor-solid cluster which 
nevertheless captures the important physics of the real system: 

~p g , r < R 

p(r) = {p s + (-p_ v -p s ) r -^,r<R (29) 
p v , r > R + w 

Xs, r < R 
X(r) = { Xs-Xs^, r<R 
0, r > R + w 

This piecewise-lincar model can be extended to give arbitrarily complex approximations and will play an important 
role below. For now, it is substituted into the expression for the free energy which is then simplified to get 

Air 

n [p] -fi v — ~y r3 i^s - W„) 

+ 4^i? 2 (wuJ (p s ,p v ,Xs) + ^~K pp (p s - p v f + -K^xs (p s - P v ) + ^-K^xl 
\ 2w w 2w , 

+ 0(Rw,w 2 ) (30) 

with 

Q v = J ui(p v ,0)dr (31) 
w (P s ,PviX) = u(p s + \ (p v - p s ) , (1 — A) x) d\ 



In the present discussion, it is assumed that the clusters are not small (in the sense that w/R < 1) so that lower 
order terms in the expression for 57 can be neglected. 

The critical cluster is a stationary point of the free energy so that the width is determined from 

= ZJ(p 8 ,p v , X s) ~ -j^ (K pp (p s - p v f + 2K™ X s (Ps - Pv) + K^xl) (32) 



and the radius from 



1 „„„ ,_ _ s2 . 1 „ px _ f - _ s , 1 M v 2 



= R K - oj v ) + 2 [wuJ (p s ,p v , Xs ) + ttK pp (p s - P v y + -K^xs (p s - P v ) + ^-K™xt (33) 
\ 2w w 2w 

= R(uj s -w„) +4wuJ(p s ,p v ,Xs) 

giving 



An, 



{P -\f } (A- (p s - ~ Pv f + 2K^ X s (P S - p v ) + K^xl) 1 (34) 



The same model applied to the vapor-liquid and liquid-solid critical clusters gives 

3V2 {<jJ v -uji) v ' 



(35) 



An. 



^ ' (Ps ' P \r (K pp (P S ~ Pf + 2K^ Xs (P S - TH) + ff~x2V 
3V2 {u)i-u s ) v ' 
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Finally, let us make the further approximation that the density of the liquid and solid states are almost the same so 
that 



64tt a; 3 / 2 fa,a,,0) 
3\/2 (oj v - Ldi) 2 
64tt uj 3 / 2 (p s ,p s ,Xs) 
3\/2 (cj; — lo s ) 2 



(K»>> (p s -p v ) 2 ) 



3/2 



(36) 



(K XX X 2 S 



3/2 



Double nucleation will be possible if oj s < oji < ui v and will be the energetically preferred pathway provided Afi„ s > 
AD, v i which in turn is more likely if one or more of these conditions is fulfilled: 

1. The barrier for a direct transition is larger than that for an indirect transitions (jj s , ~p v , Xs) > w(p ; ,p in 0) , ZU ~p vl 0) 
(It was this condition that was studied previously by Lutsko and NicolisQ). 

2. uj v — uji » u> v — u> s 

3. K PX ,K XX are not small compared to K pp . 

Clearly, the factor (uj v — loi) 2 occurring in the denominator of Ail v i can be important in raising the vapor-liquid 
barrier compared to the vapor-solid barrier. On the other hand, the factor involving the gradient coefficients is always 
going to be smaller for the vapor-liquid cluster than for the vapor-solid cluster due to the terms related to crystallinity. 

IV. GRADIENT COEFFICIENTS AND PLANAR INTERFACES 

In the standard CNT model, the excess free energy of a cluster is the sum of a bulk term and a surface term with 
the latter being proportional to the planar surface tension at coexistence. In the same way, for the model considered 
here, the planar surface tension will play a key role in determining the gradient coefficients. 



Gradient coefficients 



Since the second-gradient approximation is the result of a formal expansion of the free energy, exact expressions 
for the gradient coefficients exist (see Appendix [B]). Their evaluation requires full knowledge of the direct correlation 
function in the bulk system for all densities and crystallinitics which is of course not known. While reasonable models 
could be used to make approximate evaluations of the coefficients, the results would probably be disappointing when 
used to evaluate physical quantities such as the liquid-solid or vapor-solid planar surface tension simply because the 
squard-gradient model is a crude approximation. In fact, the relevant interfaces tend to have widths of only a few 
times the typical atomic separation so that the assumption of slowly- varying order parameters that underlies the SGA 
is unlikely to be very good, although for the particular case of the liquid-vapor interface which involves only a single 
order parameter and no structural change, it is actually rather good|l2j. Thus, while it would be interesting to make 
good evaluations of the exact expressions for the coefficients, the approach used here is more pragmatic. First, it is 
noted that a systematic expansion in crystallinity and density gives 



Kpp (p, X ) = KPP (0,0)(l + O( X ,p)) 
KP* (p,x) = C(T)p X (l + 0(x,p)) 
K xx (p lX ) = C(T)p 2 (l + 0( X ,p)) 



(37) 



with 



C(T) = 



d 2 K px (p,x) 



dpdx 



d 2 K xx (p, x) 



Second, there is evidence that the density-density coefficient can be well modeled in the liquid, i.e. for x = 0, by a 
density- independent quantity [l^. Thirdly, explicit calculations in the accessible limit of low density indicate that all 
three coefficients are relatively insensitive to the crystallinity (aside from the explicit factors shown above) at least 
up to crystallinities about half that of the bulk solid. Finally, Laird has noted[H, [36| that the structural properties 
tend to be dominated by the short-range repulsion of the pair potential so that the liquid-solid surface tension can 
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FIG. 5. Excess surface free energy for the solid-vapor, liquid-vapor and liquid-solid interfaces for a system interacting via the 
LJBG pair potential. The small symbols are simulation data taken from Broughton and Gilmer[37|] (solid-vapor and liquid- 
vapor) and from Davidchack and Laird[3(| (solid-liquid). The diamonds show the Laird approximation for the solid- liquid planar 
surface tension[3^|. The simulation results at each temperature are, from highest to lowest, for the 111, 100 and 110 planes 
respectively. 

be approximated by that of the hard-sphere solid evaluated with an effective hard-sphere diameter giving the useful 
approximation that the excess surface free energy for a planar liquid-solid interface is 

7 / s « 0.617d 2 (T) /k B T (38) 

. This suggests that K px and K xx should be dominated by hard-sphere contributions which would imply that they 
scale linearly with temperature. It is also consistent with the model free energy functional used here in which all of 
the dependence of the free energy on the crystallinity enters through the hard-sphere part of the free energy. All of 
this suggests a simple approximation for the structural coefficients along the lines of 

K? x (p : x)~C- px p X kBT<T 6 (39) 
K*x (p, X )~C xx p 2 k B Ta 6 

In a final simplification, the low-density limit gives Cp X = C xx . 

B. Planar interfaces 

Figure [5] shows the result of using this model to calculate the liquid-solid, vapor-solid and vapor-liquid surface 
tensions at various temperatures as compared to the available simulation data. It is evident that, e.g., choosing 
Cpx = Cxx t° give, e.g., a value of the liquid-solid surface tension in agreement with a single point of either simulation 
or the Davidchack-Laird approximation is enough to give a reasonable description of the liquid-solid and vapor-solid 
surface tensions over a range of temperatures. Furthermore, reasonable values for all of the physical quantities are 
found for a range of choices of the coefficient so that the model is relatively robust with respect to this choice. 

Figure [5] shows the planar profiles calculate using Cp X = C xx = 2 for a temperature below the triple point, one near 
the triple point and one above the triple point. The profiles above and below the triple point share common features: 
the width and position of the transition regions for both the average density and the crystallinity are roughly the 
same and the overall width of the interfaces is about three atomic planes. The profile near the triple point is broader 
and is actually composed of two distinct regions: the first part in which there is a modest drop in density and the 
crystallinity goes to zero, and the second region in which the density drops to that of the vapor with the crystallinity 
equal to zero. The first region has the nature of a solid-liquid transition while the second has the structure of a 
liquid-vapor transition. In fact, this can be interpreted as a caricature of a wetted interface. True wetting, with a 
liquid region of finite width, is not possible in this model because the long-ranged van der Waals forces that give 
rise to the fluid are not being explicitly modeled. Hence, the results shown are as close as this model can come to 
representing wetting - basically, wetting with the bulk fluid region having zero width. 
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FIG. 6. The order parameters (upper panels) and the 100 planar-averaged density (lower panels) for the solid-fluid interface 
at three temperatures which are, respectively, below, near and above the triple point. In the upper panels, the higher curve is 
~p(z) and the lower curve is xi z )- All quantities are in reduced units. 



V. VAPOR-CRYSTAL NUCLEATION 



A. General considerations 



Given a model for the free energy of interfacial systems it is now possible to consider the process of nucleation. 
In this context, nucleation means the formation of clusters that eventually become super-critical. As mentioned 
above, the most important issue that arises in vapor-solid nucleation is the description of the nucleation pathway, 
for which there arc at least three possibilities. The first is the conventional pathway in which the subcritical cluster 
has essentially bulk solid properties (density and crystallinity) except for an interfacial region. When the clusters 
are small, they are subcritical and the free energy increases with cluster size until a transition state - a critical 
cluster - is reached, after which further growth lowers the free energy. This pathway is therefore characterized by a 
single nucleation barrier and by near-bulk solid properties of the interior of the cluster. The second pathway involves 
double nucleation. First, a purely liquid-like cluster forms and grows until it reaches a critical size after which it is 
supercritical and stable with respect to the vapor. Within this liquid cluster, a second cluster forms, this time having 
the properties of the bulk-solid, and goes through a similar process of subcritical growth reaching a transition state 
and then becoming stable with respect to both the vapor and the liquid. This pathway is therefore characterized by 
two nucleation events and two energy barriers. The final possibility is termed transient two-step nucleation and is in 
some sense intermediate between the classical and double-nucleation scenarios. Sufficiently small clusters, which are 
of course sub-critical, are liquid-like but at a certain size, below the critical size for the liquid, the crystallinity begins 
to increase so that there is a single critical cluster, perhaps more solid-like than liquid-like, and a single energy barrier 
to be crossed on the path towards crystallization. In this case, the liquid need not even be stable with respect to the 
vapor - indeed, there is not even a priori reason why the liquid must exist as a thermodynamic state (e.g. minimum 
of the bulk free energy) at all. 

The question therefore arises as to how one uses the free energy model to characterize the nucleation pathway? 
Nucleation is of course a rare, noise-driven event and a dynamical description would seem most appropriate. In fact, in 
this sense, characterizing nucleation pathways is similar to the problem of characterizing chemical reaction pathways, 
for which the same issues occur. 



B. Double Nucleation 



The first step in characterizing any structural transition or reaction is the identification of the saddle points of 
the free energy surface, i.e. the transition states. In general, the beginning state (the bulk vapor) and the end state 
(the bulk crystal) are local minima of the free energy and the transition states are the critical nuclei which define 
the energy barrier separating the minima. Here, it is assumed that the relevant density distributions are spherically 
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FIG. 7. Phase diagram for the LJ potential with (vapor-solid) supersaturation as the independent variable. Double nucleation 
occurs in the region to the right and above the broken line: between the broken line and red line, the liquid is metastable but 
double nucleation has a higher energy barrier than single-step nucleation of the solid. 



symmetric and the order parameter fields, p (r) and \ ( r ) are again modeled by piecewise continuous functions, 

p(r) 



Pi 



Po, r < R 
+ R<r<R 
( P 2~Pi) r -^, R 



w 

w < r < R + w 



with a similar model for the crystallinity. In order to refer to different realizations of these models, the notation 
m(i,j) will be used indicating that there are i links for the density profile (which means 1 + 2i parameters since there 
is the initial radius and then one density and one width for each link) and j links for the crystallinity profile, giving 
1 + 2j parameters for a total of 2 + 21 + 2j parameters. The free energy is then a function of those parameters and the 
transition states are identified by searching the 2 + 2i + 2j parameter space for stationary points of the free energy 
surface. This is done using standard eigenvector-following techniques [2^. 

Three specific models will be studied here. First is the "CNT" model in which there is a single link in both the 
density and crystallinity profiles together with the additional constraint that the radii and widths of the two profiles 
are the same and that the density and crystallinity in the bulk region are the same as for the bulk crystal for the 
given temperature and chemical potential. This model therefore depends on only two parameters (the radius and 
width of the profiles) and is the simplest possible model. The other models studied will be m(l, 1), a single link for 
each profile, and m(2, 1), in which there are two links in the density profile. As for the planar interfaces described 
above, this is necessary to allow for wetting of the surface and will be seen to lead to a substantial reduction of the 
free energy over the single link model. On the other hand, including additional links in the crystallinity has very 
little effect. (For example, at T = 0.6 and S = 1.25, the change in free energy of the critical nucleus found using the 
to(2, 1) model and the m(2, 2) model is on the order of one percent.) 

The interesting feature of this problem is that for sufficiently high supersaturations, both the liquid and the solid are 
more stable than the vapor so that there are three transition states: one corresponding to the vapor-liquid transition, 
another to the vapor-solid transition and a third to the liquid-solid transition. All of this follows simply from the bulk 
free energies as discussed above and illustrated in FigJTJ The question addressed here is which of the possible paths 
will occur: direct transition from vapor to solid or a two-step transition via first a vapor-liquid transition and then a 
liquid-solid transition. It is assumed that whichever transition involves the lowest free energy barriers will dominate. 

Tables [J and |TT] give the relevant free energy barriers for transitions at different supersaturations for ksT = 0.5 
and 0.6 respectively. It is found that for low supersaturation, corresponding to the case of large critical nuclei, the 
barrier for the direct vapor-solid transition is lower than that for the vapor-liquid transition. However, at larger 
supersaturations, the vapor-liquid transition becomes less costly thus implying that the double-nucleation scenario is 
favored. 

Figures [S] show the structure of the critical nuclei for a few cases. In all cases, the crystallinity drops to zero before 
the density reaches that of the bulk vapor so, in some sense, the crystalline interior is separated from the bulk vapor 
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FIG. 8. The structure of the critical nucleus for T* — 0.6 and S = 1.25 as determined using the CNT, m(l.l) and m(2,l) 
parameterizations of the profiles. In each figure, the upper curve is the average density and the lower curve is the crystallinity. 
In the figure for the m(2,l) model, a dashed line is drawn at the radius where the crystallinity becomes zero: the fact that the 
second link of the density profile begins at nearly this point is a clear illustration of the role of the intermediate liquid state in 
wetting the cluster. 




FIG. 9. The same as Fig. [8] for T* — 0.5 and S = 2.08. Desipite the existence of a metastable intermediate liquid, there is no 
significant wetting behaviour in this case. 



by a liquid buffer. This is especially true if one considers that the structure is essentially liquid-like for crystallinitics 
less than about 0.3. However, at the higher temperature, a more distinctive wetting of the cluster is seen whereby 
the crystallinity drops to zero over a region in which the density drops from a solid-like to a liquid-like value followed 
by a drop of the density from liquid to vapor values. This is analogous to the planar wetting illustrated above and 
shows that this is also a feature of the critical clusters near the triple point. 



C. Transient liquid state 



Even in regions where double nucleation is not possible, the metastable liquid state could still play a role in 
nucleation. This becomes a question of the nucleation pathway which requires much more information than just the 
properties of the critical nucleus. A full description of the pathway would necessarily have to be dynamical in nature 
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TABLE I. Properties of critical nuclei for the LJBG interaction as a function of supersaturation at k B T — 0.5. The super- 
saturation and lattice density are given in the first two rows followed by the excess free energy for the nuclei for vapor-liquid 
(VL), vapor-solid (VS) and liquid-solid(LS) nucleation. The last row gives the difference in bulk free energies for the various 
transitions. Blank entries occur for cases where it was not possible to find a transition state. 
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TABLE II. Same as Table Ufor k B T = 0.6. 
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accounting for a variety of kinetic effects including heat and mass transport. It goes without saying that such a 
detailed description, while highly desirable, can be expected to be very difficult to implement. 

Alternatively, one might imagine beginning at the transition state which is a saddle point of the free energy 
functional and is characterized by a size (in numbers of atoms), N c . It would be natural then to find other stationary 
points under the constraint of fixed number of atoms in the cluster and thus to trace a path from the critical cluster 
to the bulk vapor. This is in fact similar to the methods used in simulations . However, even in the case of the 
liquid-vapor transition this gives discontinuous paths whereby the structure changes discontinuously for some value 
of N. For the liquid-vapor transition, this is a transition from a well-defined cluster with a liquid-like central density 
to a much larger structure with a central density slightly larger than the vapor. In the present case, the transition 
observed is from a well-defined crystalline cluster to one with zero crystallinity. 

An alternative, widely used in quantum chemistry and in the study of clusters, is the construction of steepest 
descent pathways away from the transition state. These are expected to be approximations to the dynamical pathway, 
especially in the case of quasi-equilibrium dynamics due to some sort of damping: an example might be the behavior 
of colloids (which can often be modeled as simple fluids) or macromolecules in solution. Many different methods are 
used to construct the steepest descent paths (also commonly called the Minimum Free Energy Pathway or MFEP) 
including heuristic methods such as the Nudged Elastic Band and the String Method. Both methods have been 
applied to nucleation problems flfl [ill but here I expand on recent work which seeks to construct the exact MFEP 
for parameterized density profiles [12|. 

All methods of constructing steepest descent paths require the notion of a distance between two points in parameter 
space since "steepest descent" literally means the path for which the energy varies most rapidly per unit distance 
in parameter space. The problem is that the various parameters used here - densities, widths, crystallinity - arc 
incommensurate. However, since they exist only to specify a density profile, a natural solution is to define a metric 
in density-space and then to use this to induce a metric in parameter space. In the study of fluids, the metric was 
taken to be the Euclidean distance between two density profiles, 

d[p 7 pf = J ' (p{v)-p'(v)fdv (40) 

It would be natural to continue to use this definition however, it is unsuitable for two reasons. The first is simply that 
it is technically difficult to evaluate. The second, more importantly, is that it leads to the metric being sensitive to 
details of the lattice structure. For example, as the radius parameter varies there is little variation in the metric until 
the radius crosses an atomic shell at which point there is very rapid variation. This is unacceptable in the present 
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context since the free energy surface constructed above is based on a separation of length scales according to which 
the order parameters vary slowly over atomic length scales. 

In order to motivate a simple alternative more in keeping with the present approach, note that for the liquid-vapor 
transition the crystallinity is identically zero, x ( r ) = Oj so that the Euclidean distance function becomes 

d[p 7 pf= f(p(r)-p'(r)fdr (41) 



which is the Euclidean distance between the K = amplitudes in the expansion of the density. In fact, the parame- 
terization of the density used here is 

p{r)=-p(v)+ P(r)x(r)e^- r + ... (42) 

jeNN 

where the first sum is over the first (nearest-neighbor) shell in wave-vector space. Hence, the quantity p (r) x ( r ) is 
the spatially-varying amplitude of the smallest non-zero wavevector. It therefore seems reasonable to treat this on a 
par with the amplitude of the zero-wavevector component and to define a distance function as 

d [p, pf = J \p (r) - ~p' (r)) 2 dr+ J (p (r) X (r) - p' (r) X ' (r)) 2 dr (43) 

which is what will be used henceforth. When the order parameters are expressed in terms of a collection of scalar 
parameters as ~p(r) = 7j(r;T) and x ( r ) = X ( r ;T) this then defines a distance between points in parameter space 

d [r, rf = j (p (r;r) - p (r ; r')) 2 dr+ J (p (r;r) x (r;r) - p (r;r') X (r;r')) 2 dr (44) 

Assuming this function is sufficiently continuous, the distance function implies a metric 

r ( dp (r ; r) dp ( r ; r) dp (r ; r) x (r ; r) dp ( r; r) x (r ; r) \ 
911 (r) = J {~bYi ar, Wi W, ) dr (45) 



The steepest descent paths are then determined by 



dr, m (0 ^ 



(46) 



ds 9 P n n , (V\ d/3n 

This equation is obviously similar to a dynamics consisting of simple relaxation driven by free energy gradients as 
discussed in Appendix [Cj However, since nucleation is, fundamentally, a noise-driven process, there is no reason to 
expect a priori that the correct path can be determined by running a deterministic dynamics backwards from the 
transition state. The idea of steepest descent is that it includes the idea of being a most probable path since it is 
in some sense the most efficient path over the barrier. (An analogy would be a multidimensional random walker in 
a potential field. In order for the walker to pass over a barrier, it must take some number of improbable steps in 
the right direction until it reaches the top of the barrier. The steepest descent path is the one involving the fewest 
number of steps. Other paths must cross the same barrier, but by including more steps, say in the "wrong" direction, 
the probability of falling backwards towards the local minimum increases.) 

The nucleation pathways have been calculated using this model for values of the supersaturation such that the 
metastable liquid is more stable than the vapor, but below the double-nucleation threshold. Figure [10] shows the 
pathway for ksT = 0.5 and S = 2.08 plotted as a function of the the total number of atoms in the cluster. It is 
clear that both the central density and central crystallinity begin at the values of the bulk vapor. This is because the 
surface tension depends on the gradient of these quantities and for very small clusters, the dominance of the surface 
tension term over the bulk free energy contribution forces the gradients to be zero. As the cluster grows, the density 
increases very rapidly while the crystallinity increases more slowly. The core of the cluster therefore densities more 
quickly than it crystallizes but the effect is minor. Figure [TT] shows the same quantities for ksT = 0.6, close to the 
triple point, and for S = 1.25. In this case, transient two-step nucleation is clearly in evidence: the density increases 
while the crystallinity remains almost at zero until the cluster consists of well over 100 atoms. These contrasting 
behaviors are compared in Fig. [T5] where the number of "crystalline" atoms, defined as 

Ncrys = [p(r)x (r) dr (47) 

Xbulk J 
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is plotted as a function of N. Near the triple point, the number of crystalline atoms does not increase appreciably 
until the cluster is over 100 atoms in size whereas it increases almost immediately at the lower temperature. This 
behavior is, incidentally, quite similar to the observations of ten Wolde and Frenkel who noted that for their model 
of globular proteins, the favored nucleation path involved two-step nucleation near the triple point but that this was 
not seen further below the triple point 0. 



VI. CONCLUSIONS 



The primary goal of this work has been they study of different pathways for the homogeneous nucleation of crystals 
from solution based on mean-field, DFT models. Simply from the bulk thermodynamics, it is possible to divide the 
phase diagram into regions for which double nucleation is and is not possible. When the bulk models are extended to 
inhomogeneous, multiphase systems - here via the squared-gradient approximation - it becomes possible to determine 
when double nucleation is more energetically favorable than single-step nucleation. Finally, by studying steepest- 
descent pathways connecting the transition states to the bulk states it was possible to illustrate transient two-step 
nucleation for the Lcnnard- Jones fluid. A summary of these investigations would be that: 
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FIG. 12. The number of crystalline atoms as a function of the total number of atoms in the cluster showing transient two-step 
nucleation near the triple point, T* = 0.6, and almost completely one-step nucleation away from the triple point. 

1. there is a lower supersaturation limit for double nucleation: at the triple point, the limit goes to one and it 
increases rapidly as the temperature is lowered. 

2. transient two-step nucleation seems to be closely tied to wetting and to therefore occur most distinctly near the 



Almost all phases of this study could be improved. The underlying DFT model is highly simplified and in particular 
the inclusion of a long-ranged van der Waals term would give a more realistic description of wetting including a finite 
wetting-layer thickness. The modeling of the gradient coefficients is crude and quite empirical. One can estimate them 
directly from the DFT model[H, [l7H19j and this might give a more realistic dependence on density and crystallinity. 
The piecewise continuous models could be used with more than the minimal number of links used here or some other, 
perhaps more physical, basis functions could be chosen (or the Euler-Lagrangc equations could be solved directly 
without parameterizing the fields). Finally, instead of the Minimum Free Energy Pathways studied here, a more 
physically meaningful approach would be to search for the most probable pathways which requires the introduction of 
dynamics and noise, but which seems quite feasible based on the approach of Heymann and Vanden-Eiinden(38j . 



I am grateful to Gregoire Nicolis for reading an early draft of this manuscript and for a number of useful suggestions. 
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The calculation of bulk thermodynamic properties for the hard-sphere solid using the given parameterization is 
complex. The reason is that the hard-sphere free energy functional diverges for p slightly greater than pi a tt - at 
typical solid densities, the divergence occurs for ~p — piatt ~ 10~ 10 . It might seem that this is no problem as one simply 
takes p = pi att however, upon closer inspection, this doesn't work. 

To explain the problem, let us consider the "correct" way to fix these two densities. In principle, an equilibrium 
state must satisfy three conditions: 



triple point. 
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Appendix A: Bulk solid properties 



dF [ Pr ] 
dF [p r ] 

dpiatt p a 

dF [ Pr ] 
9a Ptatt ,7 



= M 
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Let the solution to these equations be p* , p* u , a* . Solution of these three simultaneous equations is delicate due to 
the fact that (a) \p — pi a tt | *C 1 and (b) the surprising fact that 



dF [ Pl 



dp 



dF [ PI 



p=Pi a tfPi a tt> a 



dp 



P .Piatt. 



O(p) 



(A2) 



In words, the free energy is a very rapidly varying function of p forp near piatt- 

This situation is somewhat better than it seems since in practical calculations, we are almost always performing 
some sort of search over, or tabulation in terms of, the chemical potential. Hence, rather than being given p and 
having to solve three simultaneous equations, we can often take, e.g., piatt as the independent parameter and search 
for p* , a* satisfying 



dF [ Pr ] 



dpiatt 

dF [ Pr ] 



= 



(A3) 



P*, Piatt, a* 



da 



= 



P*, Piatt, a* 



We then calculate p given these values. (Note however that even evaluating the various derivatives numerically is 
quite delicate when the range over which ~pd 3 can be varied symmetrically might be 1CP 10 or smaller.) 

There is a heuristic which alleviates many of these problems. Intuitively, one tends to think in terms of the density 
and not to distinguish between the lattice density and the average density. In fact, in most calculations for the 
hard-sphere solid, the approximation p = piatt is made and the results appear quite reasonable. Since 



df(z,z) _ df{x,y) 



dx 



+ 



df (x,y) 



dy 



(A4) 



this suggests that it must be the case that 



dF [p r ] 



dp 



dF [p r 



P=Platt'Pla 



dp 



la it 



dF [p r 



P=Pu 



>Pi*att> a * 



dp 



(A5) 



P* >P*a 



In fact, this suspicion is borne out in practice. Combining these two points, and noting that the value of a that 
makes the free energy stationary varies slowly as a function of the various densities, an efficient practical procedure 
is to choose a value of piatt, to set p = piatt, to determine the stationary value of a and to use the left hand side of 
Eq lA5l to evaluate p. This involves a simple, controlled minimization and no need to evaluate the free energy near the 
divergence. A proof that Eq. IA5I is exact, or at least a good approximation, is to my knowledge missing and would 
be useful. 

Finally, an important point is that these technicalities play no role in the calculations presented here. Since we 
deal here with interfacial problems, the important thing is how the free energy varies as the average density varies 
from that of the solid down to a relatively low value (that of a liquid or vapor). Furthermore, in clusters, the interior 
density is always somewhat below that of the bulk solid (dwarfing the comparatively tiny difference betweenp* and 
Piatt) ■ The on ly wa y these technical points would be important would be if we tried to do a free minimization of 
the free energy for a planar interface in which case the algorithm would have to find the correct values for the bulk 
system. This would essentially mean solving Eq. IA3l which. because of the divergences very near the correct solution, 
leads to numerical challenges. However, using the piecewise-continuous approximation, we always set the bulk values 
"by hand" and thereby avoid this problem 



They arc given by 



Appendix B: Formal results for gradient coefficients 



PKlf* (r) = -± J r^cfr^r) dp ^: T)dP tv b T) d ^ 



where c(ri,r2;r) is the direct correlation function for the bulk solid. Unfortunately, the DFT models used here do 
not give realistic expressions for this quantity although some useful results are possible. In particular, expanding in 
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the crystallinity gives 

KPo =h{p)+0{ X ) (Bl) 
K<* = p X h (P) + O ( X 2 ) 

K^=p 2 h(p) + 0( X ) 

with 

h(p) = Yj Q c{ r] p)rHr (B2) 

^2 (?) = yJVi ^ }^ j c (r; p) r 3 dr 
which only requires the DCF in the bulk fluid (albeit, at solid densities). Similarly, expanding in density gives 
* a r h(r ,_J_ /■„,„ ,,„, A _ -pv(r)\ 9p(r i; r)a P (r 2 ;r) 



W (r) = y r 12 , a r 12)6 (1 - e-^j ^ ^ ^ + ... (B3) 
and combining the two expansions gives 

K"P = h(0)(l + O( X ,p)) (B4) 
K px =PXl2(0)(l + O( X ,p)) 
K^=p 2 I 2 (0)(l + O(x,pj) 

and even if they did, the coefficients evaluated from these expressions might well give poor values of the surface 
tension due to the truncation of the gradient expansion. 

Coexistence occurs at a given temperature for a single value of the chemical potential corresponding, by definition, 
to supersaturation equal to 1. 



Appendix C: Steepest Descent and Dynamics 

Typical dynamical models depend on a distinction between order parameters which are densities of conserved 
quantities (such as the total mass in the canonical ensemble) and those which are densities of non-conserved quantities 
(such as the mass in the grand canonical ensemble). We begin with the latter case of a non-conserved order parameter. 
Then, the evolution is often assumed to be of the form 

(r) _ _ T sn m 



dt SiPt (r) 

where T is a transport coefficient. If space is discretized and we denote ipt (ji) = ipu, this takes the form 

dtpti _ r <9£l [ipt] 



dt dipu 

Here, the notation indicates that Q is a function of all of the order parameters {ipu}- Let us suppose that the system 
is described by some alternate set of order parameters, {(pti}, which is complete in the sense that the two sets arc 
equivalent and the relation between them is invcrtible: <pu = <pu [ip] an d ipu = ipu [<P]- Then, it follows that 

d(t>u = r >p dfl [ipt] d(j) U . . 

dt 2f dip tj diptj 

^ dfl [<p t \ dcpti d(f>u 
depti 4-^ diptj dipt] 

Suppose we did not know about this dynamics and simply wanted to write down the steepest descent equations for 
these models. In that case, it is necessary to specify a metric. We assume that the -0-space is Euclidean so that the 
distance between two sets of fields is 

d 2 [iPt (r) , iP' t (r)] = / (ip t (r) - $ (r)) 2 dr (C2) 

2 
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This implies that 

d 2 fa (r) , 4>' t (r)] -> J] [&] - rh u K]) 2 (C3) 
Assuming sufBcient analyticity, this prescribes a Riemann geometry with metric 

^ dlp ti [<j)t}di)ti[<j)t} , n .s 
gim = ?. 777 777 ( C4 J 

In general, the steepest descent equations are 

„il d£\M 

d <Pti _ y d(j>u (C5) 

\ d<t>ta 1:1 d<j> tb J 

Comparison of this to Eq. (|C1[) shows that the dynamics is equivalent to steepest descent with the relation between 
time and the distance parameter being 

ds 

Tdt = -. — r . (C6) 

\ d<Pta 9 d<f>tb J 
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